% The main use of this function is to Compute Matrix R by wwb June 1st

function R = findR(rimS, p, Nmax, Q)
p1 = p(1); p2 = p(2); p3 = p(3);
R = zeros(8,Nmax+1);     % Should I add 1 here?
R(1,:) = 1;

for i = 2:8
    Pre_i = sum(rimS(i,:));
    for m = 1 : (Nmax+1)
        for ktilt = 1 : (Nmax+1)
            R(i,1) = R(i,1)+Q(i,m) * (1-poisscdf(m+ktilt,3)) * R(i-1,ktilt);
        end
    end 
    for k = 2 : Pre_i
        for m = 1 : (Nmax+1)
            for ktilt = 1 : (Nmax+1)
                R(i,k) = R(i,k)+Q(i,m) * poisspdf(m+ktilt-k,3) * R(i-1,ktilt);
                % I assume the parameter of poisson is 3.
            end
        end
    end
end
end